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ABSTRACT 

We study the detailed structure of galaxies at redshifts z > 2 using cosmological simulations with 
improved modeling of the interstellar medium and star formation. The simulations follow the forma- 
tion and dissociation of molecular hydrogen, and include star formation only in cold molecular gas. 
The molecular gas is more concentrated towards the center of galaxies than the atomic gas, and as 
a consequence, the resulting stellar distribution is very compact. For halos with total mass above 
10 11 Mq, the median half-mass radius of the stellar disks is 0.8 kpc at z « 3. The vertical structure of 
the molecular disk is much thinner than that of the atomic neutral gas. Relative to the non-radiative 
run, the inner regions of the dark matter halo change shape from prolate to mildly oblate and align 
with the stellar disk. However, we do not find evidence for a significant "dark disk" of dark mat- 
ter around the stellar disk. The outer halo regions retain the orientation acquired during accretion 
and mergers, and are significantly misaligned with the inner regions. The radial profile of the dark 
matter halo contracts in response to baryon dissipation, establishing an approximately isothermal 
profile throughout most of the halo. This effect can be accurately described by a modified model of 
halo contraction. The angular momentum of a fixed amount of inner dark matter is approximately 
conserved over time, while in the dissipationless case most of it is transferred outward during mergers. 
The conservation of the dark matter angular momentum provides supporting evidence for the validity 
of the halo contraction model in a hierarchical galaxy formation process. 

Subject headings: cosmology: theory — galaxies: evolution — galaxies: formation — methods: nu- 
merical — stars: formation 



1. INTRODUCTION 

With the ever increasing computer power and so- 
phistication of algorithms, dissipationless cosmological 
simulations produced a consistent picture of the large- 
scale structure formation in the ACDM cosmology (e.g. 
Springel et al. 2005; Conroy et al. 2006; Boylan-Kolchin 
et al. 2009; Tcyssier et al. 2009) as well as the detailed 
structure of individual objects (e.g. Diemand et al. 2008; 
Springel et al. 2008; Stadcl et al. 2009). In these pure 
dark matter simulations, halos have a nearly universal 
density profile with a steep inner cusp (e.g. Dubinski & 
Carlberg 1991; Navarro et al. 1996, 2010; Diemand et al. 
2005), their overall shape is triaxial (e.g. Frenk et al. 
1988; Katz 1991; Bailin & Steinmetz 2005; Allgood et al. 
2006; Bett et al. 2007), and they exhibit a complex hi- 
erarchy of substructure (e.g. Moore et al. 1999; Klypin 
et al. 1999; Ghigna et al. 2000; Diemand & Kuhlen 2008; 
Zemp et al. 2009; Vogelsberger et al. 2009; Vogelsberger 
& White 2011). 

Important missing ingredients in these dissipationless 
simulations are the physics of cosmic gas and the forma- 
tion of stars. In the large-scale structure simulations, it is 
an appropriate approximation to neglect cooling and star 
formation, since on these scales gravity is the only dy- 
namically relevant force. But on smaller scales it is nec- 
essary to include baryonic physics in order to resolve the 
internal structure of dark matter halos correctly since the 
condensation of baryons alters the phase-space structure 
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of all matter components (e.g. Dubinski 1994; Gnedin 
et al. 2004; Kazantzidis et al. 2004; Gustafsson et al. 
2006; Debattista et al. 2008; Pedrosa et al. 2010; Knebe 
et al. 2010; Duffy ct al. 2010; Tissera et al. 2010; Abadi 
et al. 2010). 

The processes of star formation and its feedback on the 
interstellar medium are complex in nature and not well 
understood. Hence, often crude, empirical prescriptions 
have been used in order to model these processes in sim- 
ulations. A common example is a Kennicutt-Schmidt- 
type recipe based on the total gas mass density pq on 
small scales of the form ps cx pQ 2 . Fine tuning is then 
achieved with the help of threshold criteria (e.g. for tem- 
perature and/or density) in order to reproduce the ob- 
served Kennicutt-Schmidt relation, Eg °t Eq 4 (Schmidt 
1959; Kennicutt 1998), where Eg is the star formation 
rate surface density and Eg is the total gas surface den- 
sity (see Schaye & Dalla Vecchia 2008 and Mayer et al. 
2008 for recent reviews). 

Recent observations indicate that the Kennicutt- 
Schmidt relation shows more complexity that cannot 
be captured by a single power law over a wide range 
of total gas surface densities (Leroy et al. 2008; Bigicl 
et al. 2008) and that the star formation rate surface den- 
sity correlates better with the surface density of molec- 
ular hydrogen, Eh 2 (Wong & Blitz 2002; Genzel et al. 
2010; Bigiel et al. 2011). The complex behavior of the 
Kennicutt-Schmidt relation as a function of the total 
gas surface density can be understood if star formation 
only happens in regions where the gas is predominantly 
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TABLE 1 



Simulation 


Non-equilibrium 
cooling 


Star 
formation 


Supernova 
metal enrichment 


Supernova 
thermal feedback 


A 


yes 


yes 


yes 


yes 


Anf 


yes 


yes 


yes 


no 


B 


no 


no 


no 


no 



in the molecular rather than atomic phase (Kravtsov 
2003; Pelupessy et al. 2006; Robertson & Kravtsov 2008; 
Pelupessy & Papadopoulos 2009; Krumholz et al. 2009; 
Gncdin et al. 2009; Gnedin & Kravtsov 2010). While this 
is still a working hypothesis, theoretical studies demon- 
strated that the transition from the atomic to molecular 
gas may also correspond to the conditions under which 
gas becomes susceptible to gravitational fragmentation 
(Krumholz et al. 2011; Glover & Clark 2011). Hence, in 
order to realistically model star formation in simulations, 
it is necessary to follow the formation and dissociation of 
molecular hydrogen in the galactic interstellar medium. 

In this paper we study the formation of high-redshift 
galaxies, using new state-of-the-art modeling of galactic 
interstellar medium and star formation. We use the phe- 
nomenological molecular hydrogen modeling introduced 
in Gncdin & Kravtsov (2011) and a star formation recipe 
based on the density of molecular hydrogen. With these 
simulations we explore the impact of baryons on the 
structure of the dark matter halos, important for many 
current observational studies. We investigate the shape 
of dark matter halos as a function of radius, important 
for gravitational lensing studies; the response of the ra- 
dial halo profile to central condensation of baryons, for 
studies of the Tully-Fisher and other galactic scaling re- 
lations; the alignment of dark halos and stellar disks; and 
the distribution and evolution of the angular momentum 
of dark matter. 

2. SIMULATIONS 
2.1. Initial conditions 

We run cosmological simulations of a periodic box with 
comoving length L\> ox = 25.6 hr 1 Mpc. We adopt a 
ACDM cosmology with a total matter density parameter 
^m,o = 0.28, dark matter density Udm,o = 0.234, baryon 
density f^B.o = 0.046, cosmological constant tt\ = 0.72, 
Hubble parameter Hq = 100/ikms -1 Mpc -1 , with h = 
0.7, linearly extrapolated normalization of the power 
spectrum erg = 0.82, and spectral index n s = 0.96, con- 
sistent with the Wilkinson Microwave Anisotropy Probe 
7- year data (Jarosik et al. 2011). 

The initial conditions were generated taking into ac- 
count a non-zero DC mode (Sirko 2005; Gncdin et al. 
2011a). The DC mode corrects for a possible deviation 
of the average matter density in the box from the univer- 
sal value, arising from the finite simulation volume. It 
is equivalent to having a constant overdensity at redshift 
z = of <5 DC ,o = Pbox,o/Puni,o - 1- In general, the DC 
mode denotes an offset of the mean of a signal or wave- 
form from zero. Usually, it is common practice to ignore 
the DC mode, but this is already a constraint on the ini- 
tial conditions and therefore the initial conditions are not 
a truly random realization. The positions and velocities 
of particles are determined by the usual Zel'dovich ap- 



proximation (Zel'dovich 1970; Klypin & Shandarin 1983; 
Efstathiou et al. 1985). 

First, we ran 5 low resolution (256 3 particles) ran- 
dom realizations of the cosmological box. In these low- 
resolution simulations, we only simulate dissipationless 
evolution with an N-body code PKDGRAV2 (Stadel 
2001). We then selected a representative box that had 
the mass function of halos at redshift z = closest to 
that expected for the whole universe (e.g., Sheth & Tor- 
men 1999). 

This box has a DC mode of Sbc.q = 0.571. In this se- 
lected box we chose to refine 7 objects in the mass range 
M 20 ob ~ 10 11 - 10 13 M Q at z = (see also Figure 2). 
A^oob is the mass within r2oob such that the enclosed 
density is 200pb, with pb being the background matter 
density. The selected objects had a quiet merger history 
after z rs 2 but apart from that, they were selected ran- 
domly by visual examination. We used the traditional 
method of refining a region of interest with a large num- 
ber of dark matter particles and leaving the rest at lower 
resolution so that we correctly account for the large-scale 
tidal forces (Katz 1991; Bcrtschingcr 2001). Around the 
7 selected objects we refined a region of 5 r2oob with an ef- 
fective dark matter resolution of 2048 3 , i.e. the high reso- 
lution dark matter particles have a mass of 1.81 x 10 5 M Q . 
In the surrounding region, we increased the dark mat- 
ter mass in buffer zones by factors of 2 3 = 8 until we 
reached the initial low resolution level. For the thickness 
of the buffer zones we chose 3 lengths of the top-level 
cells (3Lo), where L = Lb ox /256 = 143 kpc (comoving). 

The gas was initialized following the same refinement 
pattern with the baryonic power spectrum. The initial 
composition of the gas is primordial, i.e. a hydrogen 
mass fraction X = 0.76, helium mass fraction Y = 0.24, 
and metal mass fraction Z = 0. Further, we set Ann = 
1.2 x 10- 5 ^/^^/(^}B,o/^) X = 1.5 x 10~ 4 (Peebles 1993, 
equation 6-119), X H2 = 2xl0~ 6 X = 1.52xl0~ 6 (Ricotti 
et al. 2001), and X Kl = X - A Hn - Xn 2 . All the helium 
is initially in the form of Hei, i.e. Yjjci = Y, Ihcii = 0, 
F He m = 0. 

The starting redshift of the refined initial conditions 
was determined by the criterion that the root mean 
square of the density fluctuations is 0.1, which resulted 
in z lc = 108. 

2.2. Input physics 

The simulations were run with the latest version of 
the gas dynamics and A-body Adaptive Refinement 
Tree (ART) code (Kravtsov et al. 1997; Kravtsov 1999; 
Kravtsov et al. 2002; Rudd et al. 2008). The Poisson 
and fluid equations are solved in super-comoving coordi- 
nates using cell-based adaptive mesh refinement (AMR) 
techniques (Kravtsov et al. 1997, 2002). ART includes 
3-dimensional radiative transfer of ultraviolet (UV) ra- 
diation from individual stellar particles and from the ex- 
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Fig. 1. — The disk of the most massive galaxy at z pb 3 (the main halo) in simulation A in face-on (left panel) and edge-on (right panel) 
view. The molecular hydrogen is shown in red, atomic gas (H I and He i) in blue. The stars (yellow points) only form in the region where 
the molecular gas resides. The linear size of the image is « 14 kpc (physical). 



tragalactic background (modeled according to Haardt & 
Madau 2001) using the Optically Thin Variable Edding- 
ton Tensor approximation (Gnedin & Abel 2001). It in- 
cludes a non-equilibrium chemical network for hydrogen 
(Hi, Hii, and H 2 ) and helium (Hei, Hen, and He in) as 
well as non-equilibrium cooling and heating rates, which 
use the local abundances of atomic, molecular, and ionic 
species and the local UV intensity (Gnedin & Kravtsov 
2011). All these reactions are followed self-consistently in 
the course of a simulation. An empirical model for the 
formation and shielding of molecular hydrogen on the 
interstellar dust allows for more realistic star formation 
recipes based on the local density of molecular hydrogen 
(Gnedin & Kravtsov 2011). ART also includes metal en- 
richment and thermal feedback due to the Type II and 
Type la supernovae (Kravtsov 2003). 

The local star formation rate volume density ps in a 
cell is calculated as 



where pn 2 is the local mass density of molecular hydro- 
gen. The star formation time scale is given by 

T sf = min [T ff (p G ) , T raax ] ( 2 ) 

where rg = (32Gpq/3tt)~ 1 ^ 2 is the free-fall time, and 
pa the local gas density (including all hydrogen and he- 
lium species) . The maximum timescale r max is set to the 
free-fall time of gas with a hydrogen number density of 
«h = %i + ^hii + 2riH 2 = 50 cm~ 3 . The star formation 
efficiency per local free-fall time is set to eg = 0.007. 
To ensure that star formation happens only in our nu- 
merical analogs of real molecular clouds, we allow star 
formation only in cells with the molecular mass fraction 
above /h 2 = 2nn 2 /«H =0.1. These cells have a range of 
total gas density from 50 to 10 4 amu cm~ 3 for the main 
halo at z w 3 (Figure 1). Stellar particles are created 



via a Poisson process with a characteristic timescale of 
2 x 10 yr. This star formation prescription is similar to 
the recipe SF2 in Gnedin et al. (2009). 

Figure 1 shows the disc of the most massive galaxy in 
our simulation at z ks 3, which we call the main halo. 
The molecular hydrogen forms only in high density re- 
gions and hence the stars are confined to these central 
regions. In traditional star formation prescriptions based 
on the total gas density instead of the molecular hydro- 
gen density, stars would be formed over a much larger 
volume filled with the lower-density atomic gas. 

We ran three versions of the simulation, which are 
summarized in Table 1. Simulation A is a full physics 
run with radiative transfer and non-equilibrium cooling. 
Simulation Anp is the same as simulation A but with- 
out supernova thermal feedback. Metal enrichment due 
to supernovae is still included. Simulation B is a non- 
radiative version without cooling and star formation. 

In all simulations, the top level I = grid is 256 3 and 
we allow for up to 9 more refinement levels (Z max = 9) 
where each higher level is refined by a factor 2 with re- 
spect to the parent level. This results in a size of the 
smallest cells Lq = Lb ox /(256 • 2 9 ) = 279 pc (comov- 
ing). A cell is refined if its dark matter or gas mass 
exceeded 1.07 x 10 6 M Q or 1.33 x 10 5 M , respectively. 
For the dark matter, this threshold corresponds to the 
mass of about 6 high resolution particles. On each re- 
finement level I, the time step is refined as well accord- 
ing to Ai/i = Ai/o/2 1 , where Ai/q is the global time step 
on the top level mesh. The value of Av is set at the 
beginning of each top level step so that the Courant- 
Friedrichs-Lewy condition (Courant et al. 1928, 1967) is 
fulfilled on all levels (Kravtsov et al. 2002). In total our 
simulation A contains 2.89 x 10 8 dark matter particles 
and 3.89 x 10 8 gas cells at z rj 2. 

2.3. Halo selection 
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Fig. 2. — Masses of the selected halos at redshifts z ~ 4, 3 and 
2 in simulation A, as well as the final host halos at z = from the 
N-body simulation. The median mass at each epoch is indicated 
by a circle. The lower selection cutoff is a fixed mass of 10 Mq, 
which results in different number of objects at different epochs. 



For the analysis presented in this paper, we mainly 
concentrate on three snapshots at redshifts around 4, 3, 
and 2 (the exact redshifts are 3.76, 2.85, and 2.03) in 
run A. These epochs correspond to 1.69 Gyr, 2.32 Gyr, 
and 3.29 Gyr after the Big Bang in our cosmology. The 
output redshifts for different runs match to within Az = 
0.005. 

Simulation Anf at full resolution was stopped at z = 
2.77 to save computing time. We ran also a lower- 
resolution version of Anf with the 8 times more massive 
dark matter particles but otherwise the same gas physics 
and parameters. We use the z « 2 snapshot from this 
version in our analysis. 

In all snapshots we ran a variant of the Bound Density 
Maxima halo finder (Kravtsov et al. 2004) and selected 
all massive objects with M2oob > 10 11 M Q in the high 
resolution region of run A. We then matched simulation 
A with the other two simulations in order to find the cor- 
responding halos in these runs. The matching procedure 
is not a trivial task since the dynamics in the various runs 
is different and the output epochs do not match exactly 
(e.g. halos already could have merged at a given time). 
A simple matching scheme based on comoving position 
and mass is not always successful, though it gives a rea- 
sonable initial guess. We therefore check all the matches 
by eye and repair possible misidcntifications. This pro- 
cedure is rather cumbersome but absolutely necessary in 
order to get quantitatively reliable results. We only use 
halos where we could clearly identify a match in all three 
simulations. There are 12, 16, and 16 objects (from a 
total of 16, 22, and 26 objects with M 200 b > 10 n M Q ) 
at redshifts ss 4, 3, and 2, respectively, that fulfill our 
selection and quality criteria (see Figure 2). 

2.4. Spatial resolution 

The force resolution in ART is roughly 2 cell sizes 
(Kravtsov et al. 1997) but numerical relaxation processes 
could affect the structure of halos on even larger scales. 
In order to estimate the scale where numerical effects 
become important, we use the local relaxation time at 



radius r 



with 



N(r) 



tdyn(r) = 2-7T* 



GM{r) 



(3) 
(4) 



where N(r) is the enclosed number of collisionless parti- 
cles (dark matter and stars), and the dynamical time is 
set by the total enclosed matter (including gas), M(r). 
Here we use a slightly different normalization than in 
the usual expression for the local relaxation time (e.g., 
Power et al. 2003) since we dropped the factor of 8 in 
front of the logarithmic term in the denominator. Zemp 
et al. (2008) found that this normalization agrees better 
with the results of N-body simulations, where they stud- 
ied the stability of isolated high-resolution halos. Re- 
laxation processes become important on the scale r vc \ at 
time t given by the solution of t re i(r r el) = t- 

We checked for all our selected objects that r ro i « 
2 — 4 Lg. Generally, dissipative runs A and Anf have 
a smaller relaxation scale than run B, due to the higher 
number of particles within a given radius. Therefore, we 
estimate that our results are numerically converged on 
scales larger than r rcs = 4 Lg = 1.12 kpc (comoving). 
We mark the resolution scale r rcs in all the plots where 
appropriate. 

2.5. Median properties 

In the following sections, we present halo properties 
that were calculated using a profiling routine described 
in the Appendix. We calculate the median value of each 
property for all selected halos as a function of radius. 
In order to display the spread among individual halos, 
we also calculate the region between the 15th and 85th 
percentile values, i.e. the region containing 70 percent 
of the halos. This is the spread specification throughout 
the paper, unless otherwise noted. 

The median values are a useful representation because 
the selected objects are all similar, with the mass span 
of at most of an order of magnitude. The most massive 
halo in run A at the three epochs (which is not the same 
halo) has M 200 b = 5.66 x 10 11 M , 8.77 x 10 11 M Q , 
and 1.75 x 10 12 M at z w 4, 3, and 2, respectively 
(Figure 2). The median mass and the 15th and 85th 
percentiles of the 12, 16 and 16 objects at these epochs 
are M 20 ob = l-7±°i x 10 n M Q , 2.7±f 7 3 x 10 n M o , and 
respectively. The median virial radius 



3.3±|;g x 10 11 M, 

and the percentiles are r2oob = 35^ 3 kpc, 50j^5 kpc, and 



+25 



67 



7 host halos are M> 



kpc, respectively. The final masses at z = of the 



200b 



1.6±H x 10 12 M G 



3. SPHERICALLY AVERAGED PROPERTIES 

3.1. Density profile 

In order to get a first impression of the impact of 
baryon physics on the overall matter distribution, we 
show in Figure 3 the local mass density p and in Fig- 
ure 4 the local density slope 7(7*) = — dlogp/dlogr. We 
normalize the radii of each halo by ?"2oob,A, the virial 
radius of this halo in run A. With this choice we al- 
ways compare the same physical scale of a specific halo 
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Fig. 3. — Median local mass density as a function of radius for the 
dark matter (top panel) and baryons (bottom panel) at z ss 2. We 
plot the quantity pr 2 in order to reduce the dynamical range. The 
shaded areas mark the region containing 70 percent of the halos. 
The median resolution scale is marked with a vertical dashed line. 
The three simulations are described in Table 1. 

in all three different simulations and avoid complica- 
tions arising from the slightly different values of r2oob 
in the three runs. The densities are normalized by 
P2O0b = 3M2oob/47rr| oob = 200,0b- In order to reduce the 
dynamical range, we plot the quantity pr 2 in Figure 3. 

The most obvious difference is in the distribution of 
the baryons. While in the non-radiative case the baryons 
form an approximately constant density core, they form 
a steep density profile with 7 rj 3.5 in the center for the 
cases A and Anf- The central region is dominated by 
the stars. 

In the non-radiative case B, the dark matter density 
profile is roughly proportional to r _1 near the center, 
in agreement with other dissipationless simulations (e.g. 
Dubinski & Carlbcrg 1991; Navarro et al. 1996; Dicmand 
et al. 2005, 2008; Stadel et al. 2009; Navarro et al. 2010). 
In all three cases, the slope at the virial radius is 7 « 2.5. 
It is expected to be less than 3 because of the lower con- 
centration of high-redshift halos relative to the their low- 
redshift counterparts. In the runs with star formation, 
the dark matter has a roughly isothermal inner density 
profile with 7 « 2. The central dark matter density in- 
creases by over an order of magnitude compared with the 
non-radiative case. 

3.2. Enclosed mass fraction 

Figure 5 shows the local enclosed mass fraction 
M/Mtot of the different matter species in simulation A. 
The fractions sum up to unity at each radius. Gener- 
ally, the stars dominate in inner 5% of the virial radius 
(about 3 kpc), while the dark matter dominates in the 
outer parts of the halo. The gas mass near the center is 
typically higher than the dark matter mass. 

An interesting feature is the much smaller scatter of 
the combined baryon profile than of the gas or stellar 
component separately. The star formation history varies 
from object to object, but whenever the stellar density 
is higher than the median value the gas density is lower 
by a similar amount, and vice versa. This is encouraging 
for our study. While our modeling of the star formation 



Fig. 4. — Median local mass density slope as a function of radius 
for the dark matter in simulations A and B as well as for the total 
matter in simulation A, at z « 2. The shading is as in Figure 3. 

rate and feedback may not be exactly right, the effect of 
the baryon dissipation on the dark matter properties is 
calculated more robustly. 

The universal baryon fraction in our simulations is 
^b,o/^m,o = 0.164. In general, we find that the median 
baryon fraction within the virial radius, /bC^oo^a), is 
higher than universal at all times (by « 5 — 15%) for 
the runs with cooling and star formation but lower than 
universal (by « 5%) for the non-radiative case. Also, 
/B(?"200b,A) generally increases with time, up to 15% 
above the universal fraction at z = 2. The halos in 
the simulations without supernova feedback can retain 
slightly more baryons than the halos in the simulations 
with supernova feedback. But the effect is relatively weak 
(between 2% and 9%) and within the scatter among in- 
dividual halos. 

These mass fractions should only be used for a qualita- 
tive comparison between the different simulations. There 
is nothing special about our choice of r2oob,A as the 
length scale. A smooth halo profile extends much fur- 
ther (we plot it out to two virial radii), and in general 
an edge of a halo is ill defined (see also, e.g., Prada et al. 
2006; Cuesta et al. 2008). 

3.3. Concentration 

In order to describe the more concentrated matter dis- 
tribution in the simulations with cooling and star for- 
mation, we use an intrinsic and general measure for the 
concentration of a halo. It is the enclosed density within 
the radius r max (location of the peak of the circular ve- 
locity curve Vmax) in units of the critical density at z = 0: 

P{ r max) n ( Knax \ , -\ 

c v = = 2 — . (5) 

Pc,0 \-"0 r max/ 

This concentration measure has the advantage that it 
is well defined both for isolated halos and subhalos (as 
long as the peak of the circular velocity curve can be 
found) and it does not make any assumptions about a 
specific shape of the density profile (Alam et al. 2002; 
Diemand et al. 2007). In principle, cy can be derived 
from observable quantities. An alternative interpretation 
is that cy is related to the number of rotations, iV ro t, at 
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Fig. 5. — Median enclosed mass fraction Af/Aftot a * three epochs in simulation A. The median baryon fraction within r^oob A> /B(f*200b a)i 
is marked with a circle and the value is given in the plot. The median resolution scale is marked with a vertical dashed line and the shaded 
areas mark the region containing 70 percent of the halos. 




Fig. 6. — Comparison of halo concentrations cy in runs A and B Fig. 7. — Ratio of the enclosed dark matter mass in the radiative 

at the three epochs, computed using the dark matter distribution to the non-radiative runs (top panel) and the individual radial 

alone. Circles show the median concentration of all halos at a given profiles M/M (?"2oob,A) (bottom panel) at redshift z 2. 
epoch. Dashed diagonal line corresponds to the equality relation. 



r max per Hubble time, 1/H , by 

iVrot = _ Kia * = * 0.113 ^ . (6) 

Another common measure for a halo concentration is 
the virial concentration c V i r = r v i r /r s , where r s is a char- 
acteristic scale radius. The virial concentration has two 
main drawbacks: i) c v i r grows even when the inner mass 
distribution remains constant, due to the comoving def- 
inition of the virial radius, and ii) c v - lr is not well de- 
fined for subhalos. If an analytical halo density profile is 
known, it is straightforward to calculate the mapping be- 
tween cy and c V i r - In the case of an NFW profile (Navarro 
et al. 1996) see for example Figure 5 in Dicmand et al. 
(2007). However, in dissipative simulations dark matter 
no longer follows the NFW profile (the inner parts are 
modified more than the outer parts) and cy is a more 
useful measure of halo compactness. 

Figure 6 shows the concentrations of the objects in 



run A versus the concentrations in run B. Here we de- 
termine Vmax using the dark matter profile only, because 
the total matter distribution is very concentrated and 
its peak circular velocity is reached below the resolved 
scale in runs A and Anf- Gas cooling and star forma- 
tion lead to much higher halo concentrations than in the 
non-radiative case, by a factor of 10 to 10,000. We also 
find that the median concentration of the selected halos 
gradually decreases with time. 

While the concentration cy is only sensitive to the mat- 
ter distribution within r max , the processes of gas cooling 
and star formation affect the dark matter profile through- 
out the whole halo (see also Rudd et al. 2008; Abadi et al. 
2010). This is illustrated in Figure 7, where we plot the 
median enclosed mass fraction M(r)/M(r 2 oob,A)- The 
dark matter is enhanced by a factor of around 2 at 
~ 0.1r2oob,A, and by a still noticeable amount out to 
a large fraction of the virial radius. 

These results indicate that the effects of baryon con- 
densation are not confined to the regions dominated by 




Fig. 8. — Median total velocity dispersion as a function of radius 
for the dark matter (top panel) and stars (bottom panel) at z 2. 
Stars dominate the mass in the inner 0.07r2oob,Ai but the dark 
matter dispersion is affected at even larger radii. 

stellar mass and also lead to subtle (but non-negligible) 
changes in the mass distribution at larger radii. A sim- 
ilar effect is also seen in the change of shape, which we 
discuss in Section 4. 

3.4. Velocity distribution 

In Figure 8, we show the total 3D velocity dispersion 
ctot f° r the dark matter and stars. In the non-radiative 
simulation, the dark matter shows a typical velocity dis- 
persion inversion (a tot drops towards the center) that is 
well known from dissipationless simulations. However, 
in the runs with cooling and star formation the central 
velocity dispersion of dark matter increases by a factor 
4 — 5 and displays a strong negative gradient. Such a 
strong increase is due to the concentrated stellar bulge 
and disk, which dominate the gravitational potential in 
the inner regions. 

In order to gain insight into the velocity structure of 
the halos, we calculate the velocity anisotropy parameter 



with cr ra d and <7tan being the radial and tangential ve- 
locity dispersions. We find that in run B the dark mat- 
ter has the usual anisotropy profile known from N-body 
simulations (e.g. Navarro et al. 2010): close to isotropic 
in the center while becoming more radially anisotropic 
with radius, with a sharp change towards tangential 
anisotropy around r2oob,A (see Figure 9). In contrast, in 
our simulation with cooling and star formation the veloc- 
ity structure is much closer to isotropic over a large radial 
range. In the inner regions, around 0.01—0.03 T^oob.A, the 
dark matter becomes slightly tangentially anisotropic, as 
it shares the rotation of the stars and gas (and changes 
its shape, see Section 4). Only in the outer parts does 
the velocity dispersion remain radially anisotropic. 

3.5. Pseudo phase-space density 

In dissipationless N-body simulations, the combination 
of local mass density, p, and local velocity dispersion, a, 
is shown to follow a power-law over the whole resolved 



Fig. 9. — Median anisotropy parameter as a function of radius 
for the dark matter at z S3 2. 

range of radii: pa~ 3 oc r~ a , with a = 1.875 — 1.94 (e.g. 
Bertschinger 1985; Taylor & Navarro 2001; Rasia et al. 
2004; Ascasibar et al. 2004; Dehnen & McLaughlin 2005; 
Hoffman et al. 2007; Schmidt et al. 2008; Stadel et al. 
2009; Vass et al. 2009; Navarro et al. 2010). The quantity 
Q = pa~ 3 is called the "pseudo phase-space density" 
and can be generalized by using different types of the 
velocity dispersion (total, radial, or tangential) and by 
treating the power of the velocity dispersion as a free 
parameter. Schmidt et al. (2008) showed that the value 
of the slope a depends on these definitions and that, in 
general, no universal pseudo phase-space density relation 
exists. Also, the spherically averaged true phase-space 
density does not show such a perfect power law behavior 
(Stadel et al. 2009; Vass et al. 2009). 

Here we investigate how the pseudo phase-space den- 
sity of dark matter is affected by baryon dissipation. For 
illustration, Figure 10 shows the values of Q defined us- 
ing the radial velocity dispersion. In the dissipationless 
run B we confirm the power law behavior found in pre- 
vious studies. The best-fit slope between the resolution 
radius and r2oob.A is a — 1.97 at z « 2. The slope re- 
mains essentially constant with time: a = 1.96 at z rs 3, 
and a = 1.91 at z w 4, in agreement with the results of 
Hoffman et al. (2007) and Vass et al. (2009). Qualita- 
tively similar results emerge when using the tangential 
or total velocity dispersion in the definition of Q. The 
slope a progressively increases from using <7 t an to otot to 
cr ra d, similar to the findings of Hansen et al. (2006) and 
Schmidt et al. (2008). 

However, the pseudo phase-space density is dramati- 
cally reduced in the simulation with gas cooling and star 
formation, by more than a factor of 10 near the cen- 
ter. The effect is significant all the way to the virial 
radius, with the sharpest kink at s» 0.1r2oob,A- A single 
power law no longer holds in any radial range. Similar 
results were obtained also by Tissera et al. (2010). Such 
a strong reduction is due to the increased velocity dis- 
persion of dark matter in the regions dominated by stars 
and gas, despite the density also increasing (Figure 3) 
but by a smaller amount. Since the velocity dispersion 
depends on the details of star formation and feedback, 
which vary from galaxy to galaxy, our results imply that 
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Fig. 10. — Pseudo phase-space density of dark matter as a func- 
tion of radius at z 2. In order to reduce the dynamical range we 
plot the quantity pc^r 2 , similarly to Figure 3. As a guidance, the 
dashed line shows the power law with a = 15/8 = 1.875, expected 
from the secondary collapse models (Bcrtschinger 1985). 

the pseudo phase-space density profile of dark matter in 
galaxies is not universal. 

4. TRANSFORMATION OF DARK MATTER HALO 

SHAPE 

In this section, we discuss the local shape variation of 
the mass distribution, i.e. the shape as a function of 
distance r from the halo center. We use an iteration 
method (e.g. Katz 1991; Dubinski & Carlberg 1991) and 
start with a spherical shell with middle radius r. We 
take logarithmically spaced shells with 7.5 bins per dex 
in radius. Then we calculate the shape tensor S with the 
elements 

Si = Efc m fc(^-fc) a (^)j , g . 

where summation is over all particles within that shell, 
and (rfe)j is the i component of the position vector Th 
(with respect to the halo center) of a particle/cell k with 
mass mk- By diagonalizing S we obtain the eigenvec- 
tors and eigenvalues at radius r. The eigenvectors give 
the directions of the semi-principal axes. The eigenval- 
ues of S of an infinitcsimally thin ellipsoidal shell are 
equal to a 2 /3,6 2 /3, and c 2 /3, where a, 6, and c are the 
semi-principal axes, with a > b > c. Hence, the square 
roots of the eigenvalues are proportional to the lengths of 
the semi-principal axes. We then keep the length of the 
semi-major axis fixed (but the orientation can change) 
and recalculate S by summing over all particles within 
an ellipsoidal shell with semi-major axis a = r and axis 
ratios b/a and c/a, but with the new orientation. This 
iteration is repeated until convergence is reached. As a 
convergence criterion we require that the fractional dif- 
ference between two iteration steps in both axis ratios is 
smaller than 10~ 3 . When referring to the radius or dis- 
tance from the center for an ellipsoidal shape, we always 
mean the semi-major axis a. 

There are many other methods for shape determination 
used in the literature. We present tests and comparisons 
in a companion paper (Zcmp ct al. 2011), which further 
motivates our choice of procedure, since other methods 



can lead to a significant bias in the measured shape. For 
measuring the halo shape it is essential to exclude all 
particles that are in subhalos. Not removing subhalo 
particles can lead to artificially low axis ratios (spikes) 
at the location of the subhalos (see also Zemp et al. 2011). 

4.1. Axis ratios: from prolate to oblate 

In Figure 11, we show the median axis ratios b/a and 
c/a as well as the triaxiality parameter 

T ^ 1- (Wa) 2 



l-{c/af 



(9) 



for the different matter components at z w 2. Ellipsoids 
are called oblate if < T < 1/3, triaxial if 1/3 < T < 
2/3, and prolate if 2/3 < T < 1. 

In the non-radiative case B, the dark matter shows the 
well known behavior from N-body simulations where the 
halo shape is relatively round near the virial radius but 
becomes progressively prolate (i.e. c « b < a) towards 
the center, reaching b/a « 0.55 and c/a ~ 0.4 (see also, 
e.g., Bailin & Steinmetz 2005; Allgood et al. 2006; Stadel 
et al. 2009). 

In the dissipative case A, the dark matter shape is 
much rounder at the center: b/a ss 0.95 and c/a ~ 0.8, 
and the overall shape is oblate instead of prolate. Such 
transformation has been seen in previous studies (e.g., 
Katz & Gunn 1991; Evrard ct al. 1994; Dubinski 1994; 
Kazantzidis et al. 2004; Tissera et al. 2010; Lau et al. 
2011). In particular, Abadi et al. (2010) find very similar 
axis ratios in their smooth particle hydrodynamics (SPH) 
simulations: b/a « 0.95, c/a « 0.8 — 0.9. The change in 
the shape with radius is gradual and is best demonstrated 
by the triaxiality parameter T, which steadily decreases 
from T « 0.8 outside the virial radius to T s» 0.2 at the 
center. 

The baryons in run A form a strongly flattened dis- 
tribution. The gas settles into a rather thin disk, with 
c/a w 0.2 (which includes the hot and warm gas phases; 
the cold gas disk is still thinner - see Figure 1), as op- 
posed to an almost prolate elongated shape in run B 
without dissipation. The stars also form an oblate disc 
structure, with a minimum c/a &s 0.3, although it is not 
as thin as the gaseous disk. The stellar particles ex- 
perience gravitational scattering from fluctuation of the 
potential as soon as they form and therefore diffuse away 
from the plane of the disk. Near the center, the stellar 
distribution becomes less flattened as the disk transitions 
into a bulge. The shape of the stellar distribution in the 
outskirts of the halo has a large scatter because of the 
low particle number and is not meaningful. 

The shape of the total matter distribution is then a 
consequence of the combination of the individual mat- 
ter component shapes and their relative importance as 
a function of radius (shown in Figure 5). In the dis- 
sipative run A, the shape is approximately round at 
r > 0.1r2oob,A, but sharply turns oblate at the inner 
radii where baryons dominate the mass (see also Knebe 
et al. 2010). In the innermost region, the matter dis- 
tribution is determined by the more round bulge. This 
detailed shape structure of the mass distribution is im- 
portant for accurate modeling of the strong gravitational 
lensing effect due to massive elliptical galaxies, which we 
discuss in Section 7. 




Over the three epochs that we investigate in detail, we 
do not find any significant changes in the overall shape of 
the matter distribution. Generally, the shape converges 
farther from center than the density profile. For example, 
Stadel et al. (2009) found that the convergence radius for 
the shape in a pure N-body simulations was a factor 3 
larger than the convergence radius for the density pro- 
file. One should keep this in mind when interpreting the 
results in Figure 11. 

4.2. Twisting of the halo orientation 

Why does baryon dissipation change the shape of the 
dark matter distribution so drastically at the center? 
One intuitive interpretation is that dark matter particles 
respond to the flattened gravitational potential near the 
disk and transform their orbital structure. Low-angular 
momentum box orbits may be replaced by more round 
tube orbits. Indeed, we show later in Figure 15 that the 
average angular momentum of dark matter particles at 
the center increases by a factor of several relative to the 
non-dissipative run. In addition to studying individual 
particle orbits, we can conduct a test of this idea based 
on the global shape of the dark matter halo. If the oblate 
spheroid shape of dark matter follows the shape of the 
baryon disk, the orientation of the inner spheroid should 
align with that of the disk, regardless of the orientation 
of dark matter near the virial radius. 

We define the z-orientation of the baryon disk as the 



direction of the angular momentum of the gas plus stars 
in a spherical shell between 1 and 2 kpc (physical) from 
the center, JB,m, for each halo individually. We have 
chosen to exclude the central 1 kpc region in order to 
avoid possible complications due to the ill-defined angu- 
lar momentum of the bulge component. The exact radial 
range is not very important and similar ranges around 
our choice give essentially the same result. 

Figure 12 shows that the minor axis of the dark matter 
distribution is almost perfectly aligned with the angular 
momentum of the disk in the inner 4% of the virial ra- 
dius (similar alignment was seen by Colin et al. 2006). 
However the outer regions of the halo are twisted by 50° 
to 80° relative to the disk. We have confirmed that the 
orientation of the outer halo is consistent with that in 
the dissipationless run B. Therefore, it appears that near 
the virial radius the dark matter halo retains the shape 
and orientation set by the large-scale structure forma- 
tion, while in the inner regions the halo is twisted almost 
perpendicularly and aligned with the orientation of the 
baryon disk. 

This inner alignment is not just a passive response of 
dark matter to the disk formation. Roskar et al. (2010) 
show that fresh gas feeding the disk is strongly torqued 
by the hot halo gas on its way into the galaxy center. 
Since the hot gas must respond to the shape of the dark 
matter halo, the mutual alignment is a simultaneous and 
dynamic process. 
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Fig. 12. — Angle between the direction of the inner angular mo- 
mentum of the baryons Jb i n (between 1 and 2 kpc) and the direc- 
tion of the minor axis b (solid) or the cumulative specific angular 
momentum Jdm (dashed) of the dark matter at z ss 2. 
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Fig. 13. — Ratio of the enclosed mass within 2 kpc (physical) in 
runs A and B. Lines are for the main halo. For the selected halo 
sample at z ?s 4, 3, and 2 (corresponding to cosmic times 1.7 Gyr, 
2.3 Gyr, and 3.3 Gyr), we plot the median value and the error bars 
denoting the 15th and 85th percentiles. 



The alignment of inner dark matter and stars is also 
supported by the transformation of the halo velocity 
anisotropy, from radial to mildly tangential at r < 
0.1r 20 ob,A (see Figure 9). 

This flattened, rotating configuration of dark matter 
near the center resembles what Read et al. (2009) named 
the "dark disk" . In the cosmological simulations of three 
Milky Way-sized halos with the SPH code GASOLINE, 
they found a flattened dark matter structure and at- 
tributed its origin to the accretion of satellite halos along 
the directions closely aligned with the stellar disk. Our 
simulations similarly include all the cosmological accre- 
tion history, with a factor of 5 higher dark matter mass 
resolution. We do not see a clear thin disk of dark matter, 
but rather a smooth oblate distribution near the center, 
aligned with the stellar disk. Satellite halos lose most 
of their mass by tidal stripping at larger radii and do 
not contribute significantly to the inner regions near the 
disk. Read et al. (2009) also note the flattening of their 
halos, with a similar short axis axis ratio c/a « 0.8. Thus 
we agree on the shape and orientation of the inner halo, 
but our simulations do not reveal any additional "dark 
disk". The shape of the inner dark matter distribution 
is determined by the halo contraction, which we discuss 
in Section 5, and not by accretion of satellites. 

Our interpretation is also supported by the results of 
Colin et al. (2006) , who found that the growth of a galac- 
tic bar within an isolated spherical halo led to the trans- 
formation of the halo to the oblate shape corotating with 
the bar. Their simulations did not include cosmological 
satellite accretion but nevertheless found the alignment 
effect similar to our result. 

In a wider context, the alignment of the inner dark 
matter halo is important for the long-term stability of 
stellar disks. If the halo remained triaxial, it would scat- 
ter stars off the plane of the disk and thus kinematically 
heat the disk. An oblate, aligned inner halo is key to the 
survival of thin stellar disks in CDM cosmology. 

A similar picture emerges when studying the orien- 
tation of the cumulative specific angular momentum of 
dark matter, Jdm- Near the center the directions of Jdm 



and Jb,wi are well aligned, but at the virial radius they 
are separated by 47^20 degrees. Such twisting is in ex- 
cellent agreement with the AMR simulations of a sam- 
ple of about 100 galaxies by Hahn et al. (2010), who 
found the inner dark matter aligned with the stellar disk 
to within 18°, but the outer dark matter misaligned by 
ss 50° . The gradual misalignment trend can also be seen 
in the SPH simulations of Bett et al. (2010), although 
they only probed it outside 0.25 T2oob because of lower 
particle numbers. 

These results are also important for evaluating the 
accuracy of semi-analytical models of galaxy formation 
that often assume that the specific angular momentum of 
baryons equals that of dark matter and is conserved dur- 
ing the disk formation. At the virial radius the angular 
momentum vectors of the baryons and the dark matter 
are indeed approximately aligned (within 17j^g degrees), 
as was found previously in the SPH simulations of Bailin 
et al. (2005) and in the AMR simulations of Kimm et al. 
(2011). However, the memory of this orientation is com- 
pletely erased near the disk. The value of the angular 
momentum in the inner regions also changes for both 
baryons and dark matter, as we discuss in Section 6. 

5. RADIAL HALO CONTRACTION 

As we have seen so far, baryon dissipation and conden- 
sation near the center lead to a more concentrated dark 
matter distribution (e.g. Zcl'dovich et al. 1980; Barnes 
& White 1984; Blumenthal et al. 1986; Ryden & Gunn 
1987). In this section we describe this enhanced con- 
centration by the modified model of halo contraction 
(Gnedin et al. 2004, 2011b). 

First, we note that the effect of halo contraction is a 
robust outcome of our simulations, independent of any 
interpretation using our contraction model. Figure 13 
shows the increase of dark matter mass enclosed within 
a fixed physical radius (2 kpc) compared to the non- 
radiative case. The inner dark matter mass is consis- 
tently enhanced, by a factor 3 to 6. Moreover, this en- 
hancement increases steadily with time both for the main 
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Fig. 14. — Best-fit values of the parameters of the MAC model 
for all selected halos at z ss 2. Errorbars denote the 90% confidence 
limits on each parameter. Line shows the track of the main halo 
from 2Ri5.2tozft!2. The time direction is marked with an arrow. 
Dashed lines show the relation between Ao and w that gives the 
same mass enhancement factor relative to the SAC model: 100%, 
50%, and 30%. 



halo and in the median of all massive halos. In addition, 
dissipation increases the total mass within 2 kpc, baryon 
plus dark matter, by more than an order of magnitude. 

Radial halo contraction can be accurately described 
by the modified adiabatic contraction (MAC) model 
(Gncdin et al. 2004), an extension of the standard adi- 
abatic contraction (SAC) model of Blumcnthal ct al. 
(1986). The SAC model is based on conservation of 
the specific angular momentum profile of dark matter, 
j2 M ( r ) = GM{r)r. The MAC model is based instead on 
conservation of the quantity M(f)r, where f is the orbit- 
averaged radius for particles currently located at radius 



[M B>i (fi) + Mdm,^)] r t = [M BJ (f f ) + M DM ,/(r/)] r s . 

(10) 

Here Mdm,i is the initial dark matter profile, A/dm,/ 
is the final dark matter profile (and similar for the 
baryons), and rf is the final radius of the dark matter 
shell that was initially at r$. The model assumes that 
contraction of the spherically averaged halo profile can be 
described as a motion of spherical shells that do not cross 
each other. This results in Mdm,/(^/) = Mom,* (?"«)• 
Using the mass within the orbit-averaged radius f ap- 
proximately accounts for eccentricity of particle orbits in 
cosmological structure formation simulations. The orbit- 
averaged radius can be parametrized as 



r 



A 



ro 



(11) 



This power law relation reflects typical energy and ec- 
centricity distributions of particles in halos. We take the 
pivot radius to be ro = 0.03r2oob,A, which Gnedin et al. 
(2011b) showed to minimize the correlation of A and 
w. The best-fit parameters Ao and w can vary from halo 
to halo. The SAC model corresponds to Aq = 1 and 
w = 1. The amount of contraction typically increases 
with decreasing values of Aq and increasing values of w. 
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Fig. 15. — Median ratio of the cumulative specific angular mo- 
mentum squared of dark matter in runs A and B at 2 S3 2, for 
the same amount of matter. Jb is calculated at the radius r CO rr 
corresponding to the same enclosed dark matter mass as at r in 
run A, i.e. M dm ,aM = M DM ,B(7corr). 

We use a modified version of the code Contra 1 that 
determines the best-fit values for A and w by compar- 
ing the prediction to the measured dark matter profile 
in the dissipational simulation. The fitting procedure is 
described in Gncdin ct al. (2011b). 

Figure 14 shows the best-fit parameters for our selected 
halos at z 2. The majority of the halos are clustered 
around A ~ 2.3, w ~ 0.7, but some show substantial 
scatter. For each individual halo, the MAC model pro- 
vides an excellent fit of the dark matter profile, with the 
median fractional mass error of 2.6% (averaged over all 
bins). Overall, the level of contraction is weaker than 
predicted by the SAC model, confirming the conclusion 
of Gnedin et al. (2011b). 

Dashed lines in Figure 14 illustrate the strength of halo 
contraction in comparison to the SAC model. Let us de- 
fine the mass enhancement ratio at radius r: Fm{t) = 
Mr>M,f{f)/M-DM,i{f). Each line corresponds to a fixed 
fraction of this enhancement factor relative to its value in 
the SAC model: f M (r) = F M (r\A ,w)/F M (r\l,l) = 1, 
0.5, and 0.3, respectively. The radius at which these frac- 
tions are evaluated is r = 0.005 r 2 oob,A, near our resolu- 
tion limit. For most of the halos, the mass enhancement 
factor at this radius is 30% to 60% of the value expected 
in the SAC model. 

Figure 14 also shows the evolution of the main halo 
over time in the space of the model parameters. It starts 
with relatively strong contraction, reaches a peak value 
around the cosmic time 1.3 Gyr, recedes until 1.8 Gyr, 
then reaches another peak, and finally settles into a 
quasi-steady state at Ao ~ 2.4, w « 0.5. The two bouts 
of strongest contraction (when w is highest and Ao is low- 
est) do not correspond to the epoch of the peak of star 
formation following a major merger around t ~ 1.6 Gyr. 
Instead, the first bout precedes the merger and the sec- 
ond happens after the system reaches new dynamical 
equilibrium. The exact correlation with the dynamical 
state is not straightforward or very strong. At lower 
redshift the contraction effect stabilizes and is still sig- 

1 http : //www . astro . lsa. umich. edu/~ognedin/ contra/ 
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Fig. 16. — Change of the total angular momentum of dark matter 
between runs A and B at z RJ 2 for the main halo as a function of 
its enclosed mass fraction rather than radius (solid lines). Dashed 
lines show the same change for the baryons, as a function of the 
baryon mass fraction. One percent mass fraction of dark matter 
(the left boundary of the plot) corresponds to r ~ 0.007r2oob A> 
slightly larger than the resolution limit. The profile of baryons in 
run A is plotted to the resolution limit. Circles mark the values of 
the angular momentum within a fixed mass shell in run A at z RJ 2 
from Figure 17. 

nificant. 

6. EVOLUTION OF THE ANGULAR MOMENTUM 
OF DARK MATTER 

The structure of the halo is affected by the evolution of 
its angular momentum. We find that baryon dissipation 
dramatically affects the angular momentum profile of the 
inner halo. 

At a fixed radius, the total angular momentum Ldm 
(and the specific angular momentum per unit mass, Jdm) 
increases with time simply because more mass accumu- 
lates as the gas falls in and halo contracts. However, in 
the idealized scenario of particles on circular orbits and 
spherical shells that do not cross, the angular momen- 
tum of each lagrangian shell with enclosed mass My)M Is 
conserved, J dm (Mdm) = GM{r)r = const. As the en- 
closed baryon mass increases, the shell moves inward but 
retains the angular momentum profile Jum(Mdm) ■ This 
is a test of the foundation of the halo contraction model. 
In this section we consider the evolution (in time and 
between the runs A and B) of the angular momentum 
of the same spherical shells, that is, the same enclosed 
amount of material. 

Relative to the dissipationless case, the inner dark mat- 
ter shells gain a lot of angular momentum. Figure 15 
shows that the median value of Jq M (Mdm) lS greater by 
as much as a factor of 10 near the center, compared to the 
angular momentum of the same dark matter mass in run 
B. The variance between different halos is large, reflect- 
ing their different formation histories. A similar change 
of the dark matter angular momentum was also seen in 
the SPH simulation of Bett et al. (2010) at 0.1r vir , the 
innermost radius they considered. 

To further illustrate this effect, we show in Figure 16 
the change of the total angular momentum explicitly as 
a function of the enclosed mass fraction for both dark 
matter and baryons. In run A the inner 3% of baryons 
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Fig. 17. — Evolution of the cumulative angular momentum of 
the dark matter and baryons in the main halo within a radius that 
contains a constant mass. The constant mass for each species was 
chosen as the enclosed mass within 2 kpc (physical) of the halo at 
z ss 4 in run A. At z ss 2 it corresponds to 32.5% of the baryon 
mass and 2.39% of the dark matter mass (see also marked circles 
in Figure 16). 

lose most of the angular momentum (« 90%) that they 
had in run B. The loss is significant throughout the whole 
halo, out to the virial radius. In contrast, the inner 10% 
of dark matter gain angular momentum. Part of this 
gain may be due to modification of particle orbits near 
the baryon disk. Another part may be due to the low- 
Jdm particles being replaced by the high- Jdm particles 
from larger radii as a result of mergers and fluctuations of 
the potential, as first suggested by Valenzuela & Klypin 
(2003). Finally, some of the effect may be due to the 
transfer of angular momentum from the baryons to the 
dark matter. 

However, the comparison between the two runs is 
masking a strong temporal evolution in the angular mo- 
mentum profile in each simulation. A fixed amount of 
material of both baryons and dark matter, in both runs, 
systematically loses the angular momentum in the inner 
galaxy. Figure 17 shows this evolution in the main halo 
for the shell containing 2.39% of the dark matter mass 
and 32.5% of the baryon mass at z « 2 (it has radius 
of s» 1.5 kpc (dark matter) and ss 0.5 kpc (baryons), re- 
spectively, at that epoch) . The sharpest drop of Lb , by 
a factor of 2, is around the cosmic time of 1.6 Gyr. This 
epoch corresponds to a major merger event, when even 
the mass within the inner several kpc increases notice- 
ably. It causes also the largest burst of star formation 
in this galaxy, with the rate exceeding 100 M yr _1 . It 
has been known since the work of Vitvitska et al. (2002) 
that the spin of the halos decreases after mergers. A 
non-negligible fraction of particles near the center gains 
enough energy to move outside the virial radius (e.g., 
Kazantzidis et al. 2006). Later, D'Onghia & Navarro 
(2007) showed that out-of-equilibrium halos have, on av- 
erage, higher spin than relaxed systems, suggesting that 
the post-merger virialization process leads to a net de- 
crease in the halo spin. Our results show that this effect 
may operate even at radii as small as a few kpc, where 
the high-Je gas moves to larger radii and leaves the re- 
maining baryons with much lower angular momentum. 
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Fig. 18. — Angular momentum profiles of the dark matter shells 
in the main halo. For each radius r at z f» 2 we matched the same 
enclosed mass at previous epochs. Solid lines are for run A, dashed 
lines are for run B. 

Figure 17 shows that without baryon dissipation, the 
enclosed Idm(^dm) of a fixed dark matter mass Mdm 
also decreases with time. In contrast, as a result of gas 
cooling and star formation in run A, the angular mo- 
mentum of the same mass, £dm(-Mdm), remains rela- 
tively constant over time. This inner region contains a 
much larger fraction of baryon mass than dark matter 
mass, and therefore baryons carry much higher total en- 
closed angular momentum, Lb ^> £dm- The lost baryon 
angular momentum is more than enough to feed the an- 
gular momentum of the inner dark matter. Note that 
this result implies that at a fixed physical radius £dm(?") 
increases over time, because the enclosed mass Mrj^(r) 
increases as the halo contracts. 

The remarkable constancy of the angular momentum 
in run A extends to all dark matter shells out to sa 
0.1r2oot>,A- Figure 18 shows that the inner angular mo- 
mentum profile has decreased by less than 10 — 20% be- 
tween z sa 4 and z « 2. At higher z > 4, when the central 
baryon accumulation was dynamically insignificant, the 
inner £dm(-<Wdm) was fluctuating by a larger amount. 
Since z ~ 4, the dense stellar bulge and disk have sta- 
bilized the gravitational potential and helped maintain 
a very constant angular momentum profile of the dark 
matter shells. In contrast, the inner Ldm(-Wdm) of dark 
matter in the non-radiative run B decreased overall by a 
factor of several. 

Did baryons transfer just enough angular momentum 
to inner dark matter to balance the loss due to mergers? 
Details of the angular momentum evolution still need to 
be investigated in future studies, but the overall effect il- 
lustrated in Figure 18 is new and significant. The angular 
momentum profile of inner dark matter shells is fortu- 
itously constant in the region dominated by the baryons, 
at least to the final epoch of our simulation. This ef- 
fect lends strong support to the underlying equation in 



the models of halo contraction. In addition, the radial 
mixing and eccentricity of particle orbits strengthen the 
motivation for using the MAC model, which is based on 
even more robust conservation of the radial action. 

7. DISCUSSION 
7.1. Comparison with observations 

At the moment there are only a few observational con- 
straints at high redshift. van Dokkum et al. (2009) re- 
port on a massive compact galaxy at redshift z = 2.2 
with the stellar mass around 2 x 10 11 M and the line-of- 
sight velocity dispersion of 510 j^ 5 km s _1 . Converting 
this one-dimensional dispersion to the three-dimensional 
dispersion with the factor y3 gives roughly 900 km s _1 . 
Our most massive halo in run A at z sa 2 has total stellar 
mass of 1.2 x 10 11 M Q and a stellar 3D velocity dispersion 
around 725 km s" 1 (at 0.005 r 2 oob,A)- It is reasonably 
consistent with the observed analogue. 

Koopmans ct al. (2009) determine the inner slope of 
the total matter density for a sample of galaxies from the 
Sloan Lens ACS (SLACS) survey. They report an essen- 
tially constant slope 7 = 2.09^qq2 i n the inner region 
(0.2 - 1.3 R e ) in the redshift range z sa - 0.4. The me- 
dian slope for our galaxies is 7 sa 2.5 in the same region 
(Figure 4) and it does not change substantially between 
z sa 4 and z sa 2. In this inner region most of the gas has 
already been converted into stars, which dominate the 
central mass (Figure 5). If the inner density slope does 
not decrease at z < 2, as may be expected for collision- 
less systems which preserve the steepest density profile 
during a merger (e.g., Dehnen 2005; Kazantzidis et al. 
2006; Zemp et al. 2008), then our simulated galaxy pro- 
files would be too steep. However, continuous star forma- 
tion and accretion of satellites may decrease the central 
concentration and build up the extent of the stellar disk, 
thus reducing the overall density slope. The mass profile 
of the lensing spiral galaxies will soon be available from 
the new multi-band SWELLS survey (Treu et al. 2011). 

If we fit an exponential profile to each stellar disk of 
our selected objects at z sa 3 in the range between 1 
and 3 kpc, we obtain a median exponential scale length 
Rd = 470igg° P c - Since the stellar profiles are not strictly 
exponential, a more robust measure of the stellar distri- 
bution is an effective half-mass effective radius, which we 
find to be R e = 800l^Q pc. This size in reasonable agree- 
ment (and even larger) with Rafelski et al. (2011), who 
found an effective radius R e sa 600 pc for their stacked 
sample of Lyman break galaxies at z sa 3 in the same 
radial range. Our size also agrees with the Keck/OSIRIS 
integral field spectroscopy study of 12 star-forming galax- 
ies at z sa 2 — 3 by Law et al. (2009) , who found the radii 
of ionized nebula emission between 0.6 and 1.6 kpc. How- 
ever, Forster Schreiber et al. (2009) found substantially 
larger effective radii of the Ha emitting region (between 
1 and 8 kpc) in their study of 80 star-forming galaxies 
at z sa 1 - 3.5 using VLT/SINFONI integral field spec- 
troscopy. We will discuss the detailed structure of the 
disks in our galaxies in a future publication. 

Malhotra et al. (2011) showed that the half-light size 
of Lyman-a emitting galaxies does not change between 
redshift z sa 6 and z sa 2 (R e sa 1 kpc) in contrast with 
the sizes of the Lyman-break galaxies which are increas- 
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ing over the same period by a factor of 3. While both 
types of galaxy are actively star forming, the Lyman-a- 
selected objects are typically younger, relatively less mas- 
sive, and less evolved chemically. The sizes of our simu- 
lated galaxies are consistent with the Lyman-a emitters. 
Since Lyman-breaks are picking moderate-aged stars, it 
is an indication that the size of stellar disks/bulges sys- 
tematically increases with cosmic time. 

7.2. Comparison with simulations 

The problem of a very concentrated stellar distribu- 
tion is prominent in most simulations of galaxy formation 
(e.g. Govcrnato et al. 2007; Scannapicco ct al. 2009; Stin- 
son et al. 2010), although two recent SPH simulations 
seem to have overcome this by using very strong feedback 
with a supernova blastwave and applying a high star for- 
mation density threshold, for a dwarf galaxy (Governato 
et al. 2010) and a large galaxy (Guedes et al. 2011). For 
example, Governato et al. (2010) had gas cooling shut 
off in a region of 0.1 — 0.3 kpc in radius for 5 — 10 Myr 
after the supernova explosion, where most of the energy 
is lost radiatively in our simulations. 

Joung et al. (2009) performed cosmological simulations 
of galaxies at z = 3 with the AMR code Enzo. Their 
simulated galaxies are in a similar mass range to our 
sample and also show very centrally concentrated stellar 
profiles, despite substantially different treatment of gas 
physics and star formation. We have also checked that in 
our simulations the dense stellar core was already built 
up at earlier epochs, at least at z = 5. Joung et al. 
(2009) suggested possible ways around this problem: an 
early reionization to suppress gas condensation; a very 
strong feedback from stars or central black holes to re- 
duce overall star formation; or a small-scale cutoff in the 
primordial fluctuation power spectrum. 

However, Agertz et al. (2011) was able to form a more 
extended disk galaxy at z = using the AMR code 
RAMSES. They performed a parameter study with the 
star formation efficiency per free-fall time varying be- 
tween eg = 0.01 and 0.05, and found that decreasing 
eff led to a less concentrated stellar distribution, despite 
smaller stellar mass and weaker feedback. We use an 
even smaller star formation efficiency (eg = 0.007) but 
find steeper stellar profiles. This discrepancy is probably 
due to a combination of two effects. First, the Agertz 
et al. (2011) simulation has lower spatial resolution (340 
pc vs. 90 pc in our run A at z = 2) and forms stars 
based on the total gas density rather the molecular hy- 
drogen density, and at much lower density than in our 
case. Molecular hydrogen naturally forms only at high 
gas density and therefore in our simulations star forma- 
tion takes place in systematically more concentrated re- 
gions (see also Figure 1). For any given feedback model, 
the feedback is less efficient in our simulations. Second, 
Agertz et al. (2011) analyze the stellar disk at z = 0, as 
opposed to z sa 2 in our case. The high-redshift disks are 
expected to grow gradually in size and may become less 
concentrated by the present epoch. 

8. SUMMARY 

We have presented a study of the impact of baryon 
physics on the matter distribution in galaxy formation 
simulations. An important new feature of our simula- 



tions is the modeling of the physics of molecular hydro- 
gen and the star formation prescription based on the lo- 
cal molecular hydrogen density and not on the total gas 
density. We investigated the properties of the most mas- 
sive halos with Af 2 oob > 10 11 M Q to redshift z»2, using 
three simulations of the same cosmological volume that 
differ by the included gas physics and stellar feedback. 
The main conclusions of this work are: 

• The spatial distribution of the molecular gas, and 
the stars that form from it, is very compact. The 
median half-mass size of the stellar disks is 0.8 kpc 
at z ss 3. The molecular disk is much thinner than 
the distribution of the neutral atomic gas. 

• Compared to the matching dissipationless run, the 
inner regions of the dark matter halo change shape 
from prolate to oblate and align with the orienta- 
tion of the stellar disk. This alignment of the inner 
halo is important for the long-term stability of the 
thin stellar disk. 

• Despite the alignment, the inner halo is only mildly 
oblate: the short-to-long axis ratio is c/a « 0.8. 
We do not detect a "dark disk" of comparable thin- 
ness to the stellar disk. 

• The radial profile of the dark matter halo contracts 
in response to baryon dissipation. In the inner 
10% of the virial radius, the logarithmic slope of 
the dark matter density profile is constant, 7 ~ 2. 
Halo mass within the inner 2 kpc is consistently 
enhanced in all galaxies, by a factor of 3 to 6. The 
mass enhancement increases steadily with time. 
The modification of the halo profile is accurately 
described by the modified model of halo contrac- 
tion. 

• Baryon dissipation affects the structure of the dark 
matter halo far outside the extent of the stellar 
disk. The enclosed mass is enhanced by ~ 40% at 
0.2 r2oob- The halo concentration is also signifi- 
cantly increased relative to the non-radiative case. 
The halo velocity distribution is consistent with be- 
ing isotropic throughout the whole halo. 

• The pseudo phase-space density profile of dark 
matter is not universal and cannot be described by 
a single power law. In the inner region dominated 
by stars, the phase-space density is reduced by an 
order of magnitude relative to the dissipationless 
case. 

• In the non-radiative run, a fixed amount of inner 
dark matter loses most of its angular momentum 
between z — 5 and z — 2. This angular mo- 
mentum is transferred outward during major merg- 
ers and potential fluctuations. In contrast, in the 
galaxy formation run the specific angular momen- 
tum of dark matter is approximately constant in 
time, while the baryons lose most of theirs. Some 
of the baryon angular momentum is transferred to 
the inner dark matter, helping to offset the loss 
during major mergers. The approximate conserva- 
tion of the dark matter angular momentum may 
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explain why the halo contraction model describes Fermilab-KICP Supercomputing Cluster (supported by 
the modification of the halo profile so well, despite grants from Fermilab, Kavli Institute for Cosmological 
the simplicity of the model and complexity of the Physics, and the University of Chicago) and the Legato 
galaxy formation process. and Flux clusters at the University of Michigan. This re- 

search has made use of NASA's Astrophysics Data Sys- 
tem (ADS), the arXiv.org preprint server, the visualiza- 
MZ, OYG, NYC, and AVK arc supported in part by tk \ n ^ *e computer algebra system Maxima, 
NSF grant AST-0708087. The simulations and anal- and thc Python plottm S llbrar y Matplothb. 
ysis in this work have been performed on the Joint 

APPENDIX 
PROFILING ROUTINE 

In order to determine the properties of simulated galaxies, we have built a profiling routine that takes the phase- 
space coordinates from any halo finder as input. If necessary, the routine can recenter the phase-space coordinates of a 
galaxy through a shrinking sphere method (Power et al. 2003). The profiling routine then creates a radial grid between 
the specified boundaries r gr idmin and r gr id m ax- The inner boundary is set by the resolution scale of the simulation, 
and we chose a value of r gr idmin = £box/(256 • 500) = 512/500 L 9 = 286 pc (comoving). The outer boundary is set 
individually for each halo based on an estimate of the virial radius r2oob,i of halo i. Here we set r gr idmax,i = 5r2oob,i- 
For the binning, we specify a number of bins per dex in radius. This has the advantage that it is adaptive to the size 
of the individual objects: larger objects will have more bins than smaller objects. It is also more memory efficient 
than using the same number of bins for all objects. Additionally, one can choose appropriate binning depending on 
the resolution of the simulation. For the analysis presented in this paper we chose 15 bins per dex in radius. 

In order to reduce the memory load of the profiling routine, we read only parts of data at once: generally 10 7 
particles or cells. This allows us to analyze even biggest simulations on a desktop computer with a reasonable amount 
of memory. Also, the halo catalog can be split up, although it is not needed in most cases. In order to speed up the 
assignment of particles into the bins, we use a chaining mesh method and linked list data structures (e.g. Hockney &; 
Eastwood 1988). By choosing the chaining mesh size of the order of the (largest) halo size in the simulation, we only 
need to check particles in the neighboring mesh cells of a halo for being within r gr ;d ma x of that halo. 

The profiling routine calculates radial profiles for each individual matter species as well as for the total matter 
distribution. By searching for a maximum in the circular velocity curve on the grid we determine 7" max and M(r max ). 
As a quality check for a true maximum, we require that the circular velocity is smaller than the candidate maximum 
value everywhere within the extended region /check^max- Here we use /check = 1-5. 

The mass density profile of subhalos typically shows an uprise at large distances from their center due to the inclusion 
of host halo particles. We can estimate a truncation radius for subhalos by first finding where the slope 7 of the enclosed 
density profile becomes shallow enough. We take the maximum truncation radius defined by 7(r , t runcm ax) = —0.5. Using 
the enclosed average mass density has the advantage that it is much smoother than the local density profile. From 
f truncmax we then go inward in radius and find the minimum local density. The location of the minimum local density 
is used as the truncation radius r tlunc . Assuming that the background density of the host halo at the location of the 
subhalo is given by p(r trunc ), we can remove this background density and calculate the subhalo mass M(r trunc ) and 
the circular velocity curve. 
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